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I.   INTRODUCTION 

Considerable  interest  has  recently  been  shown  in  trap- 
ped waves  travelling  along  the  boundaries  of  continents.   A 
"waveguide"  effect  exists  over  the  continental  shelf.   That 
is,  wave  energy  is  confined  (essentially  by  refraction)  to 
the  continental  shelf.   Two  general  types  of  these  waves 
exist: 

A.  Edgewaves  which  are  characterized  by  wavelengths  of  hun- 
dreds of  kilometers  (km)  and  periods  of  hours  (almost  always 
less  than  a  pendulum  day) . 

B.  Shelf  waves,  which  are  generally  even  longer,  have  pe- 
riods greater  than  one  pendulum  day,  and  travel  southward 
along  the  west  coast  of  an  ocean  in  the  northern  hemisphere 
(as  do  Kelvin  waves) . 

STOKES  (1846)  showed  that  such  a  wave  guide  effect  ex- 
ists over  a  uniformly  sloping  beach  or  continental  shelf, 
with  the  amplitude  of  the  gravity  waves  decaying  exponen- 
tially to  seaward.   URSELL  (1952)  showed  that  Stoke' s  edge- 
waves  were  the  fundamental  mode  of  a  family  of  waves  order- 
ed by  the  number  of  modes  parallel  to  the  coast.   REID(1958) 
studied  long  waves  on  uniformly  sloping  shelves  of  infinite 
width,  including  the  effect  of  the  Coriolis  force.   Reid 
showed  that  the  sea  surface  may  react  as  an  "inverse  baro- 
meter" and  that  atmospheric  pressure  systems  may  be  a  driv- 
ing force  for  edgewaves.   He  found  that  the  Coriolis  force 
could  cause  the  wave  period  to  vary  from  46%  less  than  to 
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86%  greater  than  that  for  the  non-rotating  case,  depending 
on  the  direction  of  travel.   A  new  quasigeos trophic  wave  is 
now  possible,  analogous  to  a  Kelvin  wave,  having. no  small 
scale  counterpart. 

ROBINSON  (1964) initiateda  study  of  the  continental 
shelf  wave  and  studied  the  data  of  HAMON  (1962,  1963)  re- 
lating tidal  and  barometric  conditions  at  several  stations 
on  the  eastern  and  western  coasts  of  Australia.   In  this 
model  the  continental  shelf  ends  abruptly,  at  which  point 
the  depth  becomes  infinite.   He  found  that  an  inverse  baro- 
meter effect  was  exhibited  but  that  the  propagated  shelf 
waves  had  a  celerity  double  that  of  his  calculations  for 
the  western  boundary.   MOOERS  and  SMITH  (1968)  studied  the 
relation  of  sea  level  and  barometric  conditions  along  the 
Oregon  coast  for  a  period  of  nearly  one  year.  Their  statis- 
tical results  show  a  barometric  factor  of  -1.2  cm/mb  and 
predominant  sea  level  oscillations  of  0.1  and  0.35  cycles 
per  day  in  the  summer.   They  conclude  that  a  shelf  wave  of 
period  three  days  is  travelling  north.   MYSAK  and  HAMON 
(1969)  found  shelf  waves  off  the  coast  of  North  Carolina, 
but  found  no  coupling  between  the  sea  surface  and  air  pres- 
sure in  the  frequency  range  0-0.5  cpd.   ADAMS  and  BUCHWALD 
(1969)  show  that  an  equally  suitable  driving  force  for 
shelf  waves  is  the  longshore  component  of  the  geos trophic 
wind.   This  may  account  for  the  exaggerated  frequency  re- 
sponse of  the  sea  level  on  the  east  coast  observed  by 
Hamon. 
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MYSAK  (1967,  1968) extended  Robinson ' s  work,  and  discus- 
sed the  effect  of  a  continental  shelf  of  finite  width  on  the 
frequency  of  Hamon ' s  Australian  waves.   His  theoretical  so- 
lutions correspond  more  closely  to  the  observations,  al- 
though he  still  cannot  account  for  the  extremely  low  read- 
ings along  the  eastern  boundary.   He  attributes  the  discrep- 
ancy mainly  to  the  presence  of  stratified  water  and  currents 
in  the  deep  water  beyond  the  continental  shelf.   A  signifi- 
cant discrepancy  exists  between  the  dispersion  relation  and 
that  for  waves  over  an  infinitely  wide  continental  shelf. 

This  paper  is  a  study  of  the  effects  of  a  continental 
slope  and  finite  ocean  depth  upon  the  present  one-slope 
models  of  MYSAK  (1968) .   A  sharp  discontinuity  in  the  depth 
of  water  beyond  the  continental  shelf  is  not  a  common  occur- 
rence in  the  world  ocean.   It  is  interesting  to  study  the 
two-slope  situation  where  a  gently  sloping  continental  shelf 
(slope,  s  <  0.002)  and  steeper  continental  slope  (s  w  0.05) 
form  a  transition  zone  between  the  coast-line  and  deep 
water.   Three  parameters:   the  slope  of  the  continental 
shelf,  the  slope  of  the  continental  slope,  and  the  depth  of 
the  deep  water  should  have  possible  effects  on  shelf  waves. 
These  are  investigated  below. 
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II.       ANALYSIS 


h   =   s..  x(x    ^  W    ) 


TT— 


.,         T  l  V    s2x   +  W    (s    -    s2) 

,2\i  K3  (wx  *  x  ^  w2+  Wl) 


,  ,J,  ,,r  h3=  constant 

I  II  III  =    s1VI1    +    s2W2 


'— ^s  -■  \A/|   i!**!"*8  W,  *^ 


Fig.  1.   Profile  of  the  two-slope  model. 

The  model  is  characterized  by  a  gradually  sloping  con- 
tinental shelf  of  finite  width  (Region  I)  adjoining  a  steep- 
er continental  slope  (Region  II)  which  terminates  in  water 
of  uniform  depth  (Region  III) .   A  representative  slope  and 
width  of  the  continental  shelf  of   .002  and  100  km  (Mysak, 
1968a)      are  used  below.   A  representative  depth  of  the 
deep  ocean  is  5000m  and  is  the  greatest  depth  of  Region  III 
Three  slopes  will  be  used  for  the  continental  slope:  .03, 
.05  and  .08,  with  .05  used  as  the  standard  for  comparison 
with  Mysak's  model. 

The  shallow  water  equations  are  used: 

Su/at  -  f v  +  g3*yax  =  dv/at  +  fu  +  g^/3y  =  o  (1) 

a/ax(hu)  +  a/dy(hv)  +  ac/at  =  o  (2) 
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where  (u,  v)   are  the  (x,  y)   velocity  components  and   C 
is  the  free  surface  height. 

Consider  a  wave,  moving  in  the   y   direction,  speci- 
fied by 


^  =  Uk(x)ei(CTt-m^ 


vk  =  Vk(x)ei(<Jt-roy>  (3) 

Ck  =  ^(x)ei(CTt-my> 

where  k  denotes  the  region. 
Using  (1)  and  (3) ,  u  and  v  are 


ig         *~    „        _   ST?   .       i(crt-my)  ,,. 

uk   =  -2~2    (fmT7  -   cr  ^  )^e  (4) 

-g         ,__        -   ^T]   v       i(CTt-my)  ,-, 

vk   =  -2—2    (omrj  -    f  ^   )    e  (5) 

f   -cr  J 

Eliminating  u,  v  from  (2)  in  Region  I  gives  the  equa- 
tion 

2_.p2    s,  fm      ^ 

hir?i'  +  siT7i     +   (  ^T"  "     a him  )T?i  =  °  (6) 

After  making  the  substitutions 


h..  =  s,x 

2   2 
=  cr  -fz    fm 

pl     gs,  '    <J 
(6)  can  be  written 

xT^"  +  T)1'  +  (p1-m2x)T?1  =  0  (7) 
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This  equation  in  r?   has  as  its  solution  (REID,  1958) 

-z-,/2 

T^    =    lA1Fa1(z1)    +   B1Ga1(z1)}e  (8) 


where   Fa    (z    )    is   Rummer's    function,    and  Ga    (z    )    is    a    sec- 
ond   (independent)    solution    (SLATER,    1960) : 

Fa(z)    B    x   +   az    +  a(a+l)Z2    +   a(a+iya+2)z3    +    ... 

+   a(a+l) .  .. (a+n-l)zn   +    ...  (g) 

(nl)2 


and 

2 

Ga(z)    =    Fa(z)ln  z+  az  (-  -   2)+   a  (a+1)  z2  (1Ta~3f    )  ••• 

%a  a(a+l)    ' 

+   a(a+l)...(a+n-l)zn(i+^I+...  +  ^I1-2-l-...f)+... 

(10) 

m-p1 

Here,    a,    =  — ^ —   ,    z,    =   2mx,    and  A1 ,    B,    are   constants 
1  2m  1  11 

of   integration.      Since   Ga, (z1 )    approaches    infinity   as    z 

approaches    zero,    B.    must   equal    zeros 

-z,/2 
r?x    =   A1Fa1(z1)e  (11) 

Substituting    (11)    into    (4)    gives 

— zl  dFa    (z    ) 

Ul    =    -  ff^2  All<fW)Pal(zl>-20         4    1    H      (12) 

Similarly,  in  Region  II, 

-z  /2 
T?2    =   e  {A2Fa2(z2)    +   B2Ga2(z2)}  (13) 
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and 


u2  =  -^:Z2/%2[«.-f,Pa2(22,-2,^|i^] 

r  dGa  (z  )  -,, 

+  B2L(a-f)Ga2(z2)-  2cr  ^     J}         (14) 


where 


m-p2 
z2  =  2mi2    '       a2    =   IST  ' 


so  =  —   and  p„  =  -  — — 

2    s0       ^2     gs„     (7 


In  Region  III  the  counterpart  of  (6)  is 

2   2 
h3r?3"  +  (  g  gf   "  h3m2  )T73  =  °  •  <15) 


Since  T]  must  be  bounded  for  large  x, 


T]3  =  A^e     where  (16) 


2   2  -,  1 
2    cr  ~fz 


I  -  \m"   - 

L       gh 


3 


2    and 


igA?  -£(x-W  -W  ) 

u^  =  "   o   o  (fm+Cf£)e       -1  (17) 

J      cr  -f 


There  are  now  equations  defining  T)  and  U  in  each  re- 
gion.  The  next  step  is  to  patch  together  the  solutions  for 
TJ  and  U  at  the  points  x  =  W_  ,  and  x  =  (W.+W_)  thus  eliminat- 
ing the  constants  A,  ,  B,  .   The  patching  conditions  are 
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2  3 

CCJ-,  =   ^J,  =  °       (surface  height  continuity)   (18) 

r  i2  r       i3 

LhuJ,  =  LhuJ2  =  0       (normal  flux  continuity)      (19) 

where 

[    :l  ,  t    ].  -  t    ]k  . 


The  following  abbreviations  will  be  used 


F.  =   Fa  .  (z  .)     G  .  s  Ga  .  (z  .) 
3  3      3  3  D   D 


The  subscript,  1,  refers  to  the  solution  for  Region  I 
(continental  shelf)  where  it  joins  Region  II  (continental 
slope) .   The  subscript,  2,  refers  to  the  Region  II  where  it 
joins  Region  I.   The  subscript,  3,  refers  to  the  continental 
slope  where  it  joins  Region  III,  the  flat  bottom.   The  func- 
tions  subscripted  3  have  the  same  form  as  those  subscripted 
2  with  the  exception  of  the  variable,  z  ,  which  is  deter- 
mined by  the  distance  from  the  origin. 

Using  (18)  and  (19)  between  Regions  I  and  II  and  set- 
ting A  =  1,  the  constants  of  integration  A  and  B  can  be 
solved  for? 

G„F..  -  F,  G~ 

A   =  -^± ±-±-  (20) 

2    G  F  '-  F  G  ' 
2  2     2  2 

F  F  '  -  F  F  ' 
12     2  1 

B2  =      ^"T  (21) 

G2F2,_  F2G2" 
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Using  (18)  and  (19)  between  Regions  II  and  III  gives  the 
final  equation  in  terms  of  m  and  <7   only: 


G2Fl'-  F1G2' 


G2F2,_  F2G2 


-  {(4-m)F3  +  2mF3' 


F  F  '  -  F  F 
1  2     2J1 


G2F2  "  F2G2 


-  |U-m)G   +  2mG  '}  =  0  (22) 


Because  there  is  no  way  to  solve  (22)  analytically,  it 
is  necessary  to  find  the  roots  numerically  using  the  IBM 
360/67  computer  system  at  the  Naval  Postgraduate  School. 
For  a  fixed  cc  =   cr/f  a  search  routine  was  used  to  find  the 
several  m's  satisfying  the  equation.   The  computer  work  is 
described  in  Appendix  A. 
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III.   RESULTS  AND  CONCLUSIONS 

A.    REVIEW 

There  are  a  number  of  questions  to  be  answered.   MYSAK 
(1968a)  in  his  finite-width  model  founds 

1.  Shelf -wave  numbers  are  inversely  related  to  the 
shelf  width,  for  a  fixed  frequency. 

2.  There  is  a  low  wave  number  cut-off  for  edgewaves 
which  is  a  function  of  the  shelf  width.   That  is,  as  the 
shelf  width  increases,  the  smallest  possible  wave  number 
decreases  (the  largest  possible  wavelength  increases). 

3.  The  fundamental  mode  edgewaves  of  REID  (1958) 
with  periods  greater  than  one  pendulum  day  do  not  exist 
over  a  continental  shelf  of  finite  width. 

It  should  be  noted  that  Mysak's  solution  is  only  ap- 
proximate in  that  the  maximum  shelf  depth  h  is  assumed  much 
less  than  the  ocean  depth,  D,  leading  to  an  approximation 
of  the  equation  expressing  continuity  at  the  edge  of  the 
shelf  (Mysak's  equation  (10)).   His  results  (labeled  MA 
below)  depend  on  this  approximation.   Exact  solutions  of 
Mysak's  equation  (6)  (corresponding  to  equation  (22)  in  this 
work) ,  were  also  generated  so  that  comparisons  could  be 
made  among  MA,  an  exact  solution  for  Mysak's  model  (ME)  and 
results  of  the  two-slope  model  studied  in  Section  II (TS). 
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B.    COMPARISON  OF  MYSAK'S  APPROXIMATE  SOLUTION  WITH  AN 

EXACT  SOLUTION  FOR  MYSAK'S  MODEL 

Nine  cases  were  studied  in  order  to  compare  ME  and  MA. 
Two  of  the  cases  are  illustrated  in  Fig.  2.   In   all  cases, 
the  continental  shelf  is  100  km  wide  with  a  slope  of  .002, 
duplicating  Mysak's  sample  calculation.   The  bottom  depth 
varied  from  500m  to  5000m.   The  shelf  wave  results  are 
shown  in  Figs.  5-9. 


T^ > — r/i 


A    fooom 


\-^—/ookm  — »?'   )  t  } )  i ,,  / ,  /// 


Y~~-lOC  km—m^rrrXf-rrrr- 


(Vertical  exaggeration  =  100:1) 


Fig.  2.   Profile  of  the  finite-width  model 


Note  that: 

1.  The  approximate  solution  consistently  gives  a 
smaller  wave  number  for  a  particular  cc     and  mode.   It  ap- 
pears to  be  the  limiting  condition  for  the  exact  solution. 

2.  Except  for  the  fundamental  mode,  an  error  of  less 
than  1%  exists  between  corresponding  modes  of  MA  and  ME  in 
cases  where  the  depth  ratio  h/D  is  smaller  than  .067. 

3.  For  the  fundamental  mode  there  is  still  a  signifi- 
cant error  introduced  in  m  when  using  MA  for  large  depth 
ratios  and  frequencies  (>10%  for  u)  >  0.8  when  D=5000m)  . 
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4.  The  edgewave  results  are  shown  in  Figs.  18  -  20. 
Neither  MA  nor  ME  gives  a  fundamental  edgewave  mode  similar 
to  that  of  Reid.   This  can  be  seen  directly  from  Mysak's 
equation  (10)  in  the  approximate  case.   Because  the  argu- 
ment is  positive  by  definition,  a  must  be  negative  in  order 
for  the  Laguerre  function  to  have  roots  (zeroes) .   This  in 
turn  requires  that  either  f  >  a      or   a  »  f. 

5.  A  low  wave  number  cut-off  does  exist  for  ME  edge- 
waves,  which  diminishes  with  an  increase  in  depth  and  in- 
creases with  an  increase  of  |  CO  |   (Table  I)  .   The  cut-offs 
are  always  lower  than  those  for  MA,  and  are  quite  symmetric 
with  respect  to  direction  of  travel. 

Table  I.   ME  edgewave  wave  number  cut-off  (mW) 

Deep-water  depth  (m)     First  mode  Second  mode 
1000                   0.58  1.32 

3000  0.32  0.76 

5000  0.26  0.58 

As  pointed  out  by  Mysak  and  Reid,  the  edgewave  disper- 
sion relation  is  not  symmetrical  with  respect  to  direction 
of  travel,  due  to  the  influence  of  the  Coriolis  parameter. 
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C.    COMPARISON  OF  THE  MYSAK  MODEL  WITH  THE  TWO-SLOPE  MODEL 


/OOOvr, 


•  7/7? 


1     (-*-  /ookfr, 


•56  km 


/''/  7/ 


4 


300O*">n 


>  V  >  /  / 


(Vertical  exaggeration  =  100:1) 
Fig.  3.    Two-slope  model. 

The  Mysak  model  was  compared  with  a  two-slope  model 
whose  continental  slope  was  .05  (Fig.  3) .   The  results  are 
shown  in  Figs.  5-9.   Note  that: 

1.  Except  for  the  fundamental  mode  for  small  oj  ,  there 
is  little  similarity  between  the  dispersion  relations  for 
the  two  models,  especially  for  the  large  deep-water  depths 
(i.e.,  wide  continental  slopes).   For  a  deep-water  depth  of 
5000  m,  both  mode  2  and  3  waves  of  TS  have  smaller  wave 
numbers  than  mode  2  of  ME  for  any  given  frequency. 

2.  The  fundamental  TS  shelf  wave  does  not  asympto- 
tically approach  CO  =  1.0  for  large  wave  numbers  as  does  the 
corresponding  ME  wave.   In  fact,  the  fundamental  wave  now 
behaves  like  the  fundamental  edgewave  of  Reid. 

3.  A  low  wave  number  cutoff  is  still  present  for 
edgewaves,  but  at  a  significantly  lower  wave  number  (for 
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mode  1,  nearly  half  that  of  ME) .   The  cutoffs  are  no  longer 
symmetric  with  respect  to  direction  of  travel  (Table  II) . 


Table  II   Two-slope  edgewave  wave  number  cutoff  (mW..  ) 


Deep-water  depth 


CO  <  0 


co  >  0 


m 

Mode 

1 

Mode 

2 

Mode 

1 

Mode  2 

1000 

0.42 

1.32 

0.28 

1.28 

3000 

0.26 

0.76 

0.16 

0.76 

5000 

0.22 

0.58 

0.12 

0.58 

D.    COMPARISON  OF  CONTINENTAL  SLOPES 


7-ttt 


Slope  =0.03 


Slope  =0.08 


(Vertical  exaggeration  =  100  si) 

Fig.  4.   Two-slope  models  with  different  continental 
slopes. 


Three  values  were  used  for  the  value  of  the  continen- 
tal slope;  .03,  .05,  and  .08  (Fig.  4).   The  results  are 
shown  in  Figs.  10  -  18,  for  deep-water  depths  from  500  m 
to  5000  m. 
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1.  Wave  numbers  for  a  particular  mode  of  trapped 
waves  over  a  fixed  deep-water  depth  decrease  with  an  in- 
crease of  continental  slope  width. 

2.  For  the  fundamental  mode  over  a  constant  gradient 
continental  slope,  wave  numbers  increase  with  an  increase 
in  slope  width  (i.e.,  an  increase  in  deep-water  depth). 
All  other  modes  decrease  with  an  increase  in  width. 

3.  A  curve  is  not  available  for  mode  3  for  a  conti- 
nental slope  of  0.03  and  deep-water  depth  of  5000  m.   This 
is  attributed  to  unknown  problems  of  the  computer  routine. 
This  problem  does  not  occur  elsewhere. 

4.  Discontinuities  appear  in  the  dispersion  relations 
for  modes  2  and  3  with  a  continental  slope  0.05  (in  the 
deep-water  depth  range  2200-3400  m)  and  with  a  continental 
slope  0.08  (for  depths  greater  than  2800  m)  ,  suggesting  the 
presence  of  complex  values  of  m.   A  similar  phenomenon  is 
not  observed  for  a  slope  0.03.   Equation  (22)  was  investi- 
gated for  complex  values  of  m  and  real  u)  (Appendix  B)  and 
complex  roots  were  found  for  a  slope  of  0.05  and  depth 
2800  m  (Fig.  22) .   Since  the  surface  height  is  the  real 
part  of  C  (x,  y,  t)  complex  values  of  m  imply  a  spatial 
growth  rate  exp  Cpm(my) }  in  the  positive  y  direction.  The 
roots  m  are  complex  conjugates,  so  that  one  wave  grows  and 
one  decays  at  this  rate.   Then  the  most  unstable  wave 
(i.e.,  the  one  with  the  maximum  growth  rate)  would  be  ex- 
pected to  dominate  the  shelf-wave  spectrum. 
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5.  No  complex  wave  numbers  were  found  for  the  edge- 
waves  studied.  Continental  slope  width  is  inversely  pro- 
portional to  wave  number  as  in  Mysak ' s  results. 

E.    RECOMMENDED  FURTHER  STUDY 

The  next  step  would  be  to  study  further  the  effect  of 
different  continental  shelf  slopes  using  the  TS  model. 
More  study  is  mandatory  in  D  4  above,  both  to? 

1.  find  its  limits  and 

2.  determine  if  this  is  a  mathematical  curiosity  or 
a  physical  reality. 

Future  investigators  should  seek  to  avoid  approxima- 
tions to  their  models.   As  this  study  has  shown,  the  results 
for  the  exact  equations  can  be  significantly  different  than 
the  approximations. 
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Fig.  10 
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APPENDIX  A 
NUMERICAL  SOLUTION  OF  ONE-AND  TWO-SLOPE  MODELS 

//    EXEC    FORTCLGP,t>ARM.FCRT=»LIST,MAP»  ,REGI  GN.  GC=9  5K  ,T  IME.GO= 
//FORT.SYSIN    DD    * 

REAL*8    L,M,H,C0R,SIG,SL0(2)  ,W(2  )  ,ANS  ,X  (  3  ) ,  7  (  3  ) ,  P  (  2) , 
GRAV, DELTA 

COMMON  SLC, M,L,H, COR, SIG,W,X,Z, GRAV,  A,KCCL 

DIMENSION    AF(2)  ,A(3)  , ANSI (3000)  ,QM(60 ) , PT(6U ,8), KLUF(8 

DATA    PT,KLUE/4EC*G. 0,8*0/ 

READ(5,1C2)     KOMAND 

IF    KOMANO    EQUALS    CCPLANK    CARD    IN    OA^A    DECK)     THE    MYSAK 
SOLUTION    AND    APPROXIMATION    ARE    CCMPUTED.     IF    KOMAND    EQUALS 
JNY    INTEGER ( THROUGH    999)     THE    2-SLCPE    SOLUTION    IS    COMPUTED.. 

IF(KOMANO.NP.O)     GOTO    4 
CALL    MY^K 
GOTO    5CC 
A    GRAV=S.PDC2 
DELTA=C. 2000-08 
TW0PI*0ARCQS(-1.Q0>*    2. DO 
READ(5,52)    VI,COP,SLO 
WRITE(  6,56)WrCOP  ,SLO 
89    WW=W(  1)+W( 2) 
KOUNT=C 
KTPL=1 
T  Ms  fl 

H=W(  1  )*<L0(  1)    +  W(2)*<H0(?) 
CRIT=mOPI/H*0.400D-Gl 

NU^=CRIT/OELTA 

WRITE(6,62) 
KEY=0 
3    DO    2    K=1,6C 

M=DFLOAT(  KTRL)*OELT/i 

IK=1 

CS=FL0AT(K)/f.E01 

SIG=-C<*C0R 

OM(K)=CS 

DO    5    J=KTRLtNUM 

LOOK=C 

DO    1    1  =  1,2 

P( I )=( SlG*SIG-COP*COR) /(GRAV*SLC(I) )-C0P*M/$ IG 

AE(  I  )=(M-P( I ) )/(2.D0*N) 

1  A(  I  )=AF(  I) 
M 3)  =  A( 2) 

CALL    FINSIG(ANS»KEY) 
IF(KOOL.NE.C)    GOTO    ICC 
G0T0    4  1 
9    00    9    MZ=Jf5CC 
M7  1  =  M7-1 

8    ANSIIMZ  )=ANCT(M71) 
GOtq    5 

41     TF(ANS.LE.O.DO)     L0HK=1 

ANSK  J)=ANS 
21     IP(L.Le..:.<?QD-19)AN$I  (J)=U.DD 
JJ=J-1 

IE(ANSHJ)*ANSH  JJ) . LE. O.DG. AND. J.GT.l )    GOTO    30 
GO^O    5 
3C    PT(Kf  IK)  =  (ANSI(  J)/(ANSH  JJJ-ANSIU)  ) *OELT A+M ) *WW 
Ic(  IK.FO.  l.AND. J.GE.2)     KTRL=J-1 

WRITE     (6,61)    PT(K,IK) 
IF( IK.EO.IM)    GOTO    31 
IK=IK+1 
«5    M=M+DELTA 
60    FORMATdH    ,2HM=,PlC.3,2HTO,D10.3,4CX,lrHCORIOLIS/SIGMA 
*=,D10.3) 
ICC    DO    32    I¥P=IKfIM 
32    KLUF(  IMp)=K-l 

IM=IK-1 
31    WRI^E    (6,63)    DELTA, m,CS 

2  CONTINUE 
PT(6C,  1)=7.5D0G 
DO    44    LFAK=1.IM 

^4    KLUE(LEAK)=6C 
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°c 


33 
34 

35 


SCO 


00    90    KK=1,8 

MJ=KLUF(KK) 

WRITE(7,71>    H.KK.MJ 

WRITE(7,72)(PT(  II, KK)      ,11=1, MJ) 

WRITE (6,92)     VJ,(PT(II ,KK) ,II=1,VJ) 

WRITF( 6,62) 

CALL    PLr*rp(PT(l,l)  ,0M,KLUE(1)  ,1  ) 

00    33    ICK=2,8 

!F(KLUF( I SK) .LP. A)    GOTO    34 

CONTINUE 

ISK  l=fCK-l 

00    35    KRIC=2,ISK1 

C&LL    PLCTP     (PTt  1  ,KRIC)  ,QM , KLUE ( KRI C )  ,2) 

CALL    PLOTF     (PT(1,ISK    ) ,QW, KLUE ( ISK ) , 3) 

WRITE     (6,65) 

REA0(5,PP)    W(2) 

WRITE(6,56>     U,CCR,^L0 


IF( W( 2).NE.C.D0) 
STOP 


GOTP       8  9 


•52    FORM^  t(  KE1C#3) 

c3    FORMA T( 5E1C.3) 

*4    FORMAT    (  1H     ,5D2C7) 

5*    FOPVAT(  lHC,12X,lHMf19X,lHL,17X,5HSI0MA,17X,6HANSWER, 

215X,1HH,/,5D20. 7) 
56    F0PMAT(1H0,5C2G.7) 
*7    FORMAT*  1H     ,3025.12) 
*8    cqpm^  T(1HC,PX,4HC(1)  t16X,4HP(2)  ,  16  X  ,4HM  1 )  ,  1 6X ,  4HA  (  2  )  , 

416X,<HA(3)  ,/,5D20.9) 

61  P0RMAT(17H    MIN    OCCURS    AT    M=,E12.4) 

62  FOPMAT<  1H1) 

63  <=ORMAT    (18H    7EP0    OCCURS    AT    *=,F12.M 
65    FORMATS 1HC,39X,2HMW) 

71  F0RMAT(E2L.7,2I  1C) 

72  FOPMAT    (5E16.7) 
8?    FOR  WAT    (010.3) 

C2    FORMAT(  1HC,I 5,( /,5F20.7) ) 
1C2    F0PMAT(I3) 

ENO 


*URPOL'TINE    ^YSAK 

THIS    PROGRAM    SOLVFS    FOR    MYSAK'S    SOLUTION 
FDR     2-SLOPE    RFTURN    Tp    THF    MAIN    PROGRAM. 


ANO    APPROXIMATION 


89 


OATA 
RE*0 
WRIT 
WRIT 
GPA 
KOUN 
IM  =  A 
KTRL 
00  2 
CS=F 
M  =  OF 
nELT 

WP  IT 
SIG  = 
IK=1 
IKK  = 
00  5 
JJ=J 
IF(  J 
IF(  J 
LOOK 


=  1 

K=l,59 
L0*T(K)/6.E1 
LOAT( KTRL)*C.2CCO-C8 
A=C. 2000-08 
F(6,6C)    DELTA, M,CS 
CS*CCP 


J=KTRL,500 

J. EO. 500. AND. IK  .LF.IM) 
J.F0.5GO.ANO.IKK.LE. I*) 
=  0 


GOTO 

G0T0 


30 
31 
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16 
1? 


11 


13 
12 


30 
32 


31 
33 
3* 


14 


3G0 


0ELL=OELTA 

Q¥(K )=CS 

KEY=0 

A =  C. ^DC- (CQR /ST 0+( SI 0*S 

OEL=D/np 

DDELTA=   CrR*CnP*V»*W/(GP 

QMFGA=STG/C0R 
ARGU=2.DC*LAVDA 
CALL    DLAP( ANS,ARGU,£ ,NN 
CALL    ORDLAP( A ,ARGU,FPRI 
ELTA*DEL 


105 

6 


OUANT=l.cr:-KDD 
LAMDA ) 

$OT=PSOPT(QLANT 

ANCWER  =  (SQT-(L 
+  2.C0.)*DE 

4NSK  J)=4NS 

ANSJC  J)=ANSWEP 

IFdK.GT.     TV. am 

JJJ=JJ-1 

IF( J. GE. 2. AND. A 

IF( J. GE. 2. AND. A 

GOTO     5 

OI<=  =  M  +  (AN?T(  JJ) 

IF(  IK.FQ.  LAND. 

PT(K,IK)=    HTF 

IK=IK+1 

DIFl=M+(ANSJf JJ 

IF( ANSJ( JJ)*ANS 

PTT(K,IKK)=    OIF 

IKK=IKK+1 

WRITE(6,65)     DIP 

GOTO     1 

WPITE(6,64)    nic 

GO^G    1 

DIF1=M+(CNSJ( JJ 

PTT(K,IKK)=    OIF 

IKK=TKK+1 

WRITE (6,66)    DIF 
IF( IK.GT. IM.ftNp 

M=M+DELTA 

GOTO    2 

00    32    IWP=IK,    I 

KLUF(  IMP)=K-1 

IF( JJ.F0.500.AN 

GOTO     34 

00    33    IVK=IKK, 

KLU(IMK)     =K-1 

IF( IK.LT.IKK)     I 

TM=IKK-1 

CONTINUE 

00     14    KIN=1,IM 

KLU(KIN)=59 

KLUE( KIN)=59 

DO    105       w7?=ltA 

MJ  =  KLUE(M77) 

WRITE( 7f7CC)0D, 

WRTTE(7,7C1)     (P 

MJ=KLU(M77) 

WRITE(7,7C0)    DO 

WRITF(7,7C1)     (P 

00    300    KPR=1,60 

PT(KPR,V77)=cj( 

PTT(KPR.MZZ) =PT 

MJ=KLUF(MZ7) 

WRTTF(7,7CC)    DD 

WRITE (7,701)     (P 

MJ=KLU    (*Z7) 

WRITE(7t70C)    DD 

WRITE (7,701)     (P 

CONTINUE 

DO    QO       KK=1,8 


IG-CQR*CCR)/(GRAV*D/W*M) J/2.DQ 
AV*D) 

) 
) 
*(1.0D0-(0MEGA*0MFGA) )/(LAMDA* 


) 

C00/OVEGA)-OEL*(1.JDO-LODO/CVECA)  )**NS 
L*FPRT 


D. IKK. GLIM)    GOTO    2 


NSI ( JJ)* 
N$J( JJ)* 

/(AN^T  (J 
JJ.GT.2) 


) /(AN?J( 
J( JJJ) .G 

1 

,DIF1 


ANSI  (  IJJLLF. 0.00  )     OOTC    11 
ANSJ( JJJ) .LE.O.DO)     GOTO    12 

J J)-ANSI ( JJ) )*OELTA) 
KTRL=JJ-2 


JJJ)-ANr J(JJ) )*0ELT4) 
LO.  DO)  GOTO    13 


JJJ)^NSJ(JJ) )*DELTA) 


)/(ANKJ( 

1 

1 

.IKK.r.T.  I  V)    GOTO    2 


M 

P.IKK.LE. IM)     G0T0    31 

IV 

KK  =  IK 


V-7  7  ,MJ 

T     (II  ,V77)  ,11=1, VJ) 


,VZ7,MJ 
TT(II ,*7 

KPP  ,M77) 
T(KPR  ,VZ 

,VZ7tMJ 
T    (II ,V7 

,VZ7,MJ 
TKII  ,MZ 


7) ,11=1, VJ) 

*V< 
7)*W 


7)  ,11=1, vj) 
Z)  ,11=1, VJ) 
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MJ=KLUF(KK) 

M  I  I  —  K  I  I  I  I  K  K  J 

WR I  TEC  6,92)  MJJ  .  ( PTT( 1 1  , KK)  ,  I  1  =  1 ,MJ J ) 
9C  WRITF(6,92)  VJ t ( PT( I  I  ,KK)  , 1 1 =1 , MJ) 
WRITE<6,71) 

CALL  PLCTP(PTT<  1  ,1  J  ,QM,  KLU<1)  ,1 ) 
DO  35  MLK=2,7 
CALL  PLCTP(PTT( 1  ,MLK)  ,QM , KLUf KL K ) ,2 ) 

35  CONTTNUF 

CALL  °! CTP(dtt(1  ,8)  ,OVtKLU(fl)  ,3) 

WRITE(6,72) 

WRTTF(6t73) 

CALL  PLCTP(PT(1 ,1)  fOM,KLUE(ll  ,1  J 

DO    36    Mpp=2f7 

CALL    PLCTP(PT<1  ,MCP)  ,GM,KLUE(MOP),2) 

36  COMTINUF 

CALL    PLCTP(PT(1,R)  ,QM,KLUE(8)  ,3) 
WRITE(6,72) 
REAH( 5,88)    OD 
WRITC(f,56)    Q,OD,COR,W 

if(do.ne.c.oc)   com  *9 

RF^UPN 

c2    FORMAT (4E10.  31 

5  3    FORMAT! 5E 10. 3) 

54    FORMAT    ( 1H     ,5D2C.7) 

5C    FORMAT I  1HC,12X,1HM,1<5X,1HL,17X,5HSIGMA,17X,6HANSWER,  15 

*502L.7) 

56  FOPMftT(  1Hg,402C7) 

57  priPMAT(  1H     ,302?. 12) 

5  8    FORMA  T(  1HC,8X,^-HP(1)  ,16X,4HP(2),16X,4HA(1),16X,4HA(2), 

4/,502C9) 
60    F0pm&T(3H    M=»D1C.3,2HT0,0U.3  ,-OX , 15HS I GMA/COR IOL IS=, 

0 1 C  .  3 ) 
62    FOR*AT(  1H1) 

64  POPMAT(35H    THE     APPROXIMATION    HAS    A    ROOT     AT    M=,F15.6,4() 

H.THE    EQUATION    COES    NOT    HAVE    A    ROOT    HERE.) 

65  P0PMAT(35H    THE    APPROXIMATION   HAS    A    ROOT    AT    M=,E15.6,27 

H.THE    EQUATION    HAS    A     ROOT    AT,E15.6) 

66  C0RMAT(27H    THE    EQUATION    HAS    A    ROOT    AT , E 15 . 6 ,49H.THE    AP 

PROXIMATICN    DOES    MOT    HAVE    A    SOLUTICN    HERE.) 

71  FORMAT* 1H1,27X,26HPL0T    OF    MYSAK'S    ENTIRE    EQN) 

72  FORMA T(lHCt4CXf2HMW) 

73  FORMAT( 1H1,25X,29HPL0T    OF    MYSAK'S    APPROXIMATION) 
SB    FORMAT(01C3) 

=2    FORMA T(lHCtI5t(/.5E20. 7)1 
7CC    ^0PMAT(F2C.7,2I 1C) 
7C1    FORMATf 5E12.4) 

ENO 

SUBROUTINE    FINSIG(AN*,KEY) 

REAL* 8  L,M,X(3) ,P(3) ,FPRI (3) ,G(3),GDRI (3) ,  C, H,GR AV, COR 

,SIG,W<2),  SL0(2) ,ANS.Z<3) 

COMMON  SLO,M,L,H,COP,SIG,W,X,Z,GRAV,  A,KOCL 

OTMENSICN  AC  3) 

k  n  n  i  =  p 
0( M*M*GPAV*H+COP*COP-SIG*SIG)/ (GRAV*H) 
IF  (C  .LF.C.OOGOTO  11 
L=DSORT(C ) 
X(  1)  =  W( 1) 

X(2)=SLC(  1)*M1) /SL0(2) 
X(3)=X(2)+W(2) 
00    2    1=1,3 
It  I  )=Xf  I)*2.0C*M 
CALL    OLAP(F( I)  ,7(1 )  ,A(I)  .NN) 
CALL    DLADGOfFU  )  ,  7  ( I  )  ,G  (  I )  .  A  (  I )  ) 
CALL    DROLAP(A(I  )  ,Z(I)  ,FPRI  (I)  ) 
CALL    DGDLAP(A(I  )  ,Z(I  )  ,FDRI  (I)  ,G(I)  ,GDRnn  ) 

IF(F( I ).GE.0.103G.OR.G(I).GE.0.1030.0R.Frpi (I ).0E. 
50.103C.OR.GPRIU  )  .GE.C.103  J)    GHTO   4 
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CONTINUE 

*N*=(L-M)*(F(3)*(F(1)*GPPI  (2)-G(2)*FPRI(l)H-G(3)*(F(2) 
2*FPRId)-Ed)*FPRI<2)  )H-2*M*(FPRI  (3) *(F(1)*GPRI (21- 
3G(2)*FPRIf II J+GPPI (3)*(F(I)*FPRI (l)-F(l)*FPPI(2) )) 


1 
11 


51 
61 


WRITE(6,51) 
100*)    GOTO    4 


TF(KEY.EO.l) 

IFUNS.GE.'J 

RPTIjPN 

L=9.99n-2C 
GHTfi     1 
*OOL=i 
WRITE (6,611 

GHTR      1 


FnpMATdw     ,3D2C.?I 
FnpMAT(3«HGSCLUTI0N 


ANS 


HGFS  NCT  EXIST  REYOND  M=,E12.4I 


END 


SUBROUTINE    OLAPf Y,X,£ ,NN) 

DL*F  SOLVES  FOP  LA  GUERRE  FUNCTIONS  Of1 
ALLY  THF  ARGUMENT,/ f WOULD  RE  EXPECTED 
IS    POSITIVE, NN    IS    SET    TO    1. 

RFAL*«  Y,X,YY,YX,Y7,YK,YN 

YN  =  A 

YY=1.0C 

YK=1.DC 

YX=1.00 

YZ=1.0C 

I F  (  A  )  3,1,2 

2  MN=1 

3  YX=X/YZ*YN/YZ*YX 
YY=YY+YX 
YN=YN+1.D0 


THE  CIPST  KIND.GENER- 
TO  BE  NEGATIVE. WHEN  A 


Y7  =  Y7d.DG 

IF(DABS( YXI.GT.0.50E-C7)    GOTO 
1     Y  =  YY 
RETURN 

5C    F0RMAK1H     ,  3  (  5X  ,E  14.  7  )  | 
51    FORMAT    ( 1H     ,5020.7) 


END 

SLBROLTINE    OLA PGG ( Y, X ,DL , A ) 

DLAPGG    SOLVES    FO"    LAGUEPRE    FUNCTIONS    OF    THE    SECOND    KIND. WHEN 
A     APPROACHES    A     NEGATIVE     T MTEGER , CCNTR CL    IS    TRANSFERRED    TO 
DLAPG.WHEN    CONVERGENCE    DOE*1    NOT    ^AKE    °LACE    BEFORE    AN    ANSWER 
IS    OBTAIN-0, CONTROL    IS    RETURNED    TO    THE    CALLING    ROUTINE    AND    A 
MF^SAGC    PRINTFD;WHEN     X=C  ,     THIS    FUNCTION    I*     I NDET CR MINATE. 

PE,AL*fl    DN,DQ,DX,D2,X,Y,DL,DA 
KOUNT=0 

IF(X.EO.C)    GCTO    6 
DN=1.00 
D0=1.DC 
DX=l.DC 
DL=Y*DLOG( X) 
DA=A 
D2=0.DC 
4    DX  =  DX*DA  /DN*X/DN 

D2=D2  +  (  1.DC/DAI-I2.DC/DN) 

D0=DX*D2 

DL=DL+DO 

DA  =  DAd.DC 

1/DA     APPROACHING    INFINITY? 
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IF(DABMDA>.LE.G.1D-G2>    GOTO    1 

ANSWER    TOO    LAPGF    WITHOUT    CONVERGENCE? 

IF    (n4BS(DX).GE.l.D65.0R.OABS(nX).LE.1.0-65 >    GOTO    1C 

CONVFRGENCE? 

IF(DABS(DQ).LE.C.5D-C8. AND.K0UNT.GT.2)    GOTO    3 
DN=DN+1.00 
KOUNT  =  KOUNT-H 
GOTH    A 

1  CALL    DLAPG(Y,X,DL,A) 
GOTO    3 

2  WRITE* 6,54>KGUNT 

3  RETURN 

1C    WRITE(6,57)     KOUNT,DX,0L 
GOTO    2 
6    WPITF(6,53) 
GOTO     3 

c3    FQRMATfcSH    X    IS    C.L^GLERRE    FUNCTION    DF    THE    SECONH    KIND 
DOFS    NCT    EXI  ST) 

54    FOPMAT(  1H     ,1C0X,I3) 

*S    F0RMAT(1H     ,602C7, /f02C.7| 

57    F0RMAT(2CH       OVF°FLOW    APPROACHING. AFTER, 13 « 15W    ITER- 
ATIONS ,DX=  ,E 12 . 5 ,9H    ANO    PL  =  ,E12.5) 

END 


SUBROUTINE    DLAPG(Y,X,OL,A) 

OLAPG    SOLVES    FOR    LAGUERRE    FUNCTIONS    OF    THE    SFCCNO    KINC    WHEN 
A     IS    A    NEGATIVE     INTEGFP.Tn    AVOID    DIVIDING    BY     ZERO, THE    NTH 
TFPM    IS    THF     SUV    OF    N-l    MULTIPLICATIONS. 


22 


21 

56 


PEAL*S 

KOUNT= 

0N=1.D 

D0=1.D 

DX=1.D 

DL=Y*D 

DA  =  A 

0Z=1.D 

02=0.0 

dx=ox* 

011  =  C 
DO    1     I 

01=1. 
K0OK=I 
DD2  K  = 
IF(K.E 
IF(DAB 
01=01* 
IF(K00 
D^  3  L 
IF(DAB 
D1=D1* 
011=01 

D2=D2 
OQ=DQ* 
IF(0O. 
7Z=(D1 
KOUNT= 
IF(KOU 
DL=DL+ 
IF(DAB 
ON=ON+ 
GHTO  2 
RETURN 
WPITEC 
G0T0    2 


Y,X,D4,Z7,DZ,D1 ,02,PN,0L,Pll,DXf DO 

0 

c 
c 

0 
LOG( X) 

0 

c 

x 

.DC 

=C,KOUNT 

00 

♦  1 

CI 

Q.I)    GOTO    4 

S(D1).LE.C.50-1C)     GOTH    i 

(OA+DFLOATCK) ) 

K.GT.KCUMT)     GnTC    1 

=KOCK,KOUNT 

S(D1).LE.0.5D-1C)    GOTO    1 

(D4+DrL0AT(L) ) 

1+01 

-( 2. 00 /ON) 

DN*DN 

GE.C. 1D40)    GOTO    56 
1    +D1*(DA+DFL0AT(K0UNT)  )*02)*DX/DQ 
KOUNT+1 

NT.GT.25)    GOTO    56 
ZZ 
S(ZZ) .LE.C.50-C8)    GOTO    21 

l.DC 
2 

6,53)    ZZ 
1 
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51    FQPMATllH     ,M5X,E15.8)) 

53    FOPMAT(36H    G(X)    DID    NCT    CONVERGE.    LAST    TFRM    WAS,E15.8> 

END 

SUBROUTINE    DRDLAP(A,X, RESULT) 

DRDLAP    SOLVES    THE    FIRST    DERIVATIVE    OF    THE    LAGUEPRE    FUNCTION 
OF    THE    FIRST    KIND. 

REAL*?    YDA,YPD, RESULT, YDB,YDN,X 

YDA*A*1.DG 

RESULT=A 

YDN-l.DO 

YDP=A 

22  YOB=YDN+l.Dj 
YDD=YDD*YDA*X/    (YDN*YDB) 
RFSULT=PF^ULT+YDD 

IHDABM  YCD).LT.0.5D-8)     GOTO    21 
YDA=YDA-H.DO 
YDN= YDN+l.nO 
GOTO    2  2 
21    RETUPN 

23  FORMAT( 1H    ,2(5XtE15.Q) ) 
54,    FORMAT    (  1H    ,6020.7) 

END 

SUBROUTINE    DGDLAP(  A  ,  X, RESULT  ,  PLG  ,RESIIL  ) 

DGDLAP    SOLVES    THE    FI&ST    DERIVATIVE    OF    THE    SECOND    KIND    OF 
LAGUERPE    FUNCTION. WHEN    X=0,    THIS    FUNCTION    DOES    NOT     EXIST. 

PEAL*8     X,RESLLtRESULT,XA,PLG,XN,X3  ,X4  ,X2,XX,X?,XC 

IF(X.EO.C)    GCTO    1C 

XN=2.DG 

XA=A 

RESUL=RESULT*    OLOG ( X) +PLG/X    +1 . CO- { 2 . DO *X A ) 

X2=1.DG/XA-2.DC 
XX=XA 

21  XA  =  X^+1.DC 
IF(DABS(XA).LE.C.5D-M    GCTO    22 
XX=XX*X/( XN*XN-XN) *X* 

X2=X2    -2.DC/XN    +1.D0/XA 

XZ=XX*X2 

RESUL  =  PEHL    +  XZ 

IF(DABS( X7).LE. 0.50-9)    GOTO    10 

XN=XN    +1.DC 

IF    (DABS( X7) .GE.1.D65)    GCTO    23 

GO    TO    21 

22  KISS=1 

XN=2.DC 

XA=A 

RESUL=PE$ULT*DLOG< X)    +■    RLG/X    +  1  .  00-(  2.  DC*X  A  ) 

X0=1.DC 

XX=1.DC 

X2=-2.DC 
5    XX=XX*X 

X3=0.DC 

DO    6    11=0, KISS 

XA-l.DO 

KIN=KISS+1 

DO    7    JJ=C,II 

IF(JJ.EO.II)    GOTO    9 

IF(OABS( X4).LE.G. 50-10)    GOTO    6 
7    X4=XA*(XA+DFL0AT(  JJ)  ) 
9     IF(KIN.GT.KISS)     GOTO    6 

DO    8    LL=KIN,KISS 
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IF(PABS( X4).LE.0.50-1C)    GOTO    6 
8    X4=X4*(  XA+OFLOAT(LD) 
6    X3=X3+X4 

X2=X2-(2.PC/XN) 

XO=XO*XN'*DFLCAT(KI  SS) 

XZ=( X3+X4*( XA+O FLOAT ( KTSS) ) *X2) *XX/XQ 

KISS=KISS+1 

RFSUL=PESUL+X7 

IFCPABSI X7).LE.C.5P-8)    GOTO    10 


XN=XN+1.D0 
IF(KISS.LT.25) 
WRITE(6,54>     XZ 

23    N=XN 

WRITE(6,24)     N 

10    RFTUPN 


GOTO 


24    FORM^T(ACH    OVERFLOW    ABOUT    TO    OCCUR     IN    OGOLAP    AFTeRtI3, 

2    7H    TERMS.) 
51    FOPMATf  1H     f4(5X.E15.8) ) 
54    FORMA  T(1H     ,36HGMX)    PIP    NOT    C CNVFRGF. LAST    TERM    WAS, 

1F15.B) 

ENO 
//GO.FT06FGG1  PO  OCB= ( PECFM  =  F* ,B LKS I  7 E=l 3  3 ) ,S PACE= ( CYL , ( 15, 1 
//GP.SYSIN  00  * 
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APPENDIX  B 

NUMERICAL  SOLUTION  FOR  TWO-SLOPE  COMPLEX  ROOTS 

//  EXEC  FORTCLGP, PARM. FORT= • LI  ST, MAP* , REGION. GO= 100K, T IME . GO 
//FORT.SYSIN  DD  * 

COMPLEX* 16  M,L,ANS,X(3),Z(3),P(2),A(3),AE(2),ANSI (1620 
REAL*8  H,CnR,SIG, SLO( 2) , W( 2 ) , GRAV, PARI, PAR 2, ANSWER ( 162 
COMMON  M,L,X,Z,A,H,COR,SIG,SLO,W, GRAV,KOOL 
DATA  ANSI/1620*(0.E0,0.E0)/ 

THIS  PROGRAM  SOLVES  FOR  THE  COMPLEX  ROOTS  OF  THE  TWO 
SLOPE  MODEL.  IN  THIS  CASE, THE  ROOT  WAS   FOUND  FOR  OM- 
EGA=12/60  AT  J=  10  AND  K=  30. THIS  GIVES  AN  ANSWER  OF 
(  1.5459r-07,0.1249E-07)  FOR  THE  ROOT. 

GRAV=9.8D2 
READ(5,52)  W,COR,SLO 
WRITE(6,56)W,C0R, SLO 
39  WW=W( 1)+W(2) 

H=W(1)*SL0( 1)  +W(2)*SL0(2) 

WRITE(6,100) 

PARl=1.541234E-07 

SIG=-12.00*C0R/6.00E0  1 

DC  2  J=1,2C 

PAR1  =  PAR1  +  0. 5000-10 

M=PAR1  *(1.00,0.D0) 

DO  1  1=1,2 

P(  I)  =  ( SIG*SIG-COR*COR) /(GRAV* SLO ( I) )-COR*M/SIG 

AE(I)=(M-P( I )  )/(2.D0*M) 

1  A(I)=AE( I) 
A(3)=A(2) 

CALL    F1NSIG(ANS,KEY) 

ANSI ( J)=ANS 

ANSWER(J)=CnABS(ANSI( J) ) 

DC    2    K=221,260 

PAR2=    DFLOAT(K)*    C.50CD-10 

M=PAR1*(1.CO,O.DO)    ♦    PAR2*    (O.DO,l.DO) 

DO    8    1=1,2 

P(I)=( SIG*SIG-CQR*COR)/(GRAV*SLO( I) )-CQR*M/SIG 

AE(I)  =  (M-P( I )  )/(2.D0*M) 

8  AU)  =  AE(  I) 
A(3)=A(2) 

CALL    FINSIG(ANS,KEY) 
KCFE=(K-22C)*20    +    J 

PRINTS    ANSK21)-    ANSK820) 

ANSI(KOFE)=ANS 

ANSWER (K0FE)=CDA8S(ANS I (KOFE) ) 

M=M-(2.DC*PAR2*(0.D0, l.DO) ) 

DO  9  1=1,2 

P( I)=( SIG*SIG-CCR*COR)/(GRAV*SLC( I) )-COR*M/SIG 

AE(  I  )=(M-P( I  )  )/(2.D0*M) 

9  A(I)=AE( I) 
A(3)=A(2) 

CALL  FINSIG(ANS,KEY) 
KOFF=(K-180)*20  +  J 

PRINTS  ANSM821)-  ANSI!  1620) 

ANSI (KOFF)=ANS 

ANSWER ( KOFF)=CDABS (AN SI (KOFF) ) 

2  CONTINUE 

DC  4  LPA=1,31 
K1=(LPA-1)*20+1 
K2=LPA*20 

WRITE (6, 101) (ANSI (LP) , ANSWER (LP)  ,LP=Ki,K2) 
4  CONTINUE 

52    FORMAK5E10.3) 

56    FORMAT(lH0,5D2C7) 

100  FORMATdHl  ) 

101  FORMAT( 1H0,(9D14.5) ) 

STCP 
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END 

SUBROUTINE  F INS IG ( ANS ,KEY ) 

COMPLEX* 16  M,L,X(3),Z(3),A<3)  ,F ( 3 ) ,FPRI ( 3 ) ,G(3),GPRI 
*  (3), CANS 

RFAL*8  H,COR,SIG,SLO(2),W(2) ,GRAV, CRAPS 

CCMMON  M,L,X,Z,A,H,COR,SIG,SLO,W,GRAV,KOOL 
KOCL=0 

CRAPS=0. 1D30 

C=(M*M*GRAV*H+COR*CGR-SIG*SIG)/(GRAV*H) 

L=CDS0RT(C) 

X(  1)=W( I) 

X(2)=SLO( 1)*W( l)/SLO( 2) 

X(3)=X(2)+W(2> 

DO  2  1=1,3 

Z(I)=X(I )*2.DO*M 

CALL  DLAP(F( I  )  ,Z(  I)  ,A< I) 9NN) 

CALL  DLAPGGtF ( I ) ,Z( I ) ,G( I ),A( I ) ) 

CALL  DRDLAP( A(I ) , ZCI) ,FPRI< [J ) 

CALL  DGDLAP(A(I ) ,Z(I)  ,FPRI( I )  ,G( I ) ,GPRI< I ) ) 

IF(CDABS(FU  )  )  . GE. CRAPS. OR. CDABS(GU )) .GE .CRAPS. OR    CD- 

ABS(FPRI ( I )) .GE.CRAPS.OR.CDABSCGPRK I ) ).GE. CRAPS)    GOTO 
?    J 

2  CONTINUE 

ANS=(L-M)*{F(3)*(  F(1)*GPRI(2)-G(2)*FPRI<1)  )*G(3)*(F<2) 
2*FPRI ( l)-Fll)*FRRI(2)  ) ) +2*M* ( FPRI  ( 3 ) * ( F( 1 ) *GPRI( 2 I- 
3G(2)*FPRI (1) )+GPRI<3)*(F< I ) *FPR I <  1 )-F ( 1)*FPRI<2) ) ) 
IF(KEY.EQ.l)  WRITE(6,51)  ANS,M,L 
IF(CDABSIANS) .GE. O.IDA)  GOTO  4 
1  RETURN 
4  KOGL=l 

WRITE(6,61)    M 
GCTO    1 

51    FORMATdH     ,3D20.7) 

61  FCRMAT(3AH0SGLUTIGN  DOES  NOT  EXIST  BEYOND  M=,E12.4) 

END 

SUBROUTINE  OLAP( Y  ,  X, A,NN) 

CLAP  SOLVES  FOR  LAGUEPRE  FUNCTIONS  OF  THE  FIRST  KIND. GENER- 
ALLY THE  ARGUMENT, A, WOULD  BE  EXPECTED  TO  BE  NEGAT I VE .WHEN  A 
IS  POSITIVE, NN  IS  SET  TO  1. 

C0MPLEX*16  Y,X,A,YN,YX 

REAL*8  YY,YZ 

YN=A 

YY=1.D0 

YK=1.DC 

YX=1.D0 

YZ=1.C0 

3  YX=X/YZ*YN/YZ*YX 
YY=YY+YX 
YN=YN+1.D0 
YZ=YZ*1.D0 

IF<COABS(YX) .GT.0.5D-07)  GOTO  3 
1  Y  =  YY 
RETURN 

50  FCRMATdH  ,  3(  5X, E 14. 7  )  ) 

51  FORMAT  (1H  ,5020.7) 

END 


SLBROUTINE  DL APGG ( Y, X ,0L , A) 
CCMPLEX*16  Y, X,DL,A,DA,DQ,DX,D; 
REAL*8  ON 
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CLAPGG  SOLVES  FOR  LAGUERRE  FUNCTIONS  OF  THE  SECOND  KINO. WHEN 
A  APPROACHES  A  NEGATIVE  INTEGER, CONTROL  IS  TRANSFERRED  TO 
DLAPG.WHEN  CONVERGENCE  DOES  NOT  TAKE  PLACE  BEFORE  AN  ANSWER 
IS  OBTAINED, CONTROL  IS  RETURNED  TO  THE  CALLING  ROUTINE  AND  A 
MESSAGE  PRINTED;WHEN  X=0,  THIS  FUNCTION  IS  INDETERMINATE. 

KOUNT=0 

IKCDABS(X).EQ.O.DO)  GOTO  6 
DN=l.DO 
DQ=l.CO 
DX=l.DO 
DL=Y*CDLOG(X) 
DA=A 
D2=O.DO 
<i    OX=DX*DA/DN*X/DN 

D2  =  D2+(l.DO/DA)-(2.r>0/DN) 

D0=DX*D2 

DL=DL+DU 

DA=DA+l.DO 

1/DA  APPROACHING  INFINITY? 

IF(CDABStCX)  .GE.0.lD65.OR.CDABS(DX  )  .L E .0 . 1D-60 )  GOTO  1 
ANSWER  TOO  LARGE  WITHOUT  CONVERGENCE? 

IF(CDABS(DA) .LE.0.1D-2)  GCTO  1 

CCNVERGENCE? 

1F(CDABS(DQ) .LE. 0 . 5D- 8. AND. KOLNT. GT .2  )  GOTC  3 

DN=DN+1.D0 

K0UNT=K0UNT+1 

GOTO  A 

1  CALL  CLAPG(Y,X,DL,A) 
GOTO  3 

2  WRITE(6,54)K0UNT 

3  RETURN 

10  WPITE(6,57)  KCUNT,DX,CL 
GOTO  2 
6  WPITE(6,53) 
GCTO  3 

53  FCRMAT(59H  X  IS  0. LAGUERRE  FUNCTION  OF  THE  SECOND  KIND 
H    CCES  NOT  FXIST) 

54  FCRMATUH     ,100X,I3) 

55  F0RMAT(1H    , 6D20 . 7 , /, D20 . 7  ) 

57  FCRMATI29H   OVERFLOW  APPROACH ING . AFTER  ,I3,15H  ITERAT- 
/  I0NS,DX=,E12.5,9H  AND  DL=,E12.5) 

END 

SUBROUTINE  DL APG( Y, X, CL , A ) 
COMPLEX* 16  Y,X,A,DL,DA,ZZ,D1,D11,DX 
REAL*8  D2,DN,DQ 

DLAPG  SOLVES  FOR  LAGUERRE  FUNCTIONS  OF  THE  SECONC  KIND  WHEN 
A  IS  A  NEGATIVE  INTEGER. TO  AVOID  DIVIDING  BY  ZERO, THE  NTH 
TERM  IS  THE  SUM  OF  N-l  MULTIPLICATIONS. 

KOUNT=0 
DN=1.D0 
DQ=1.D0 
DX=1.D0 
DL=Y*CDLOG<X) 
DA=A 
DZ=1.C0 
D2=0.D0 
22    DX=DX*X 

C11=0.D0 
DO  1  I=0,KCUNT 
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D1=1.D0 
KC0K=I+1 
D02    K=0,I 
IF(K.EQ.I  )    GUTO    4 
IF(CDABS(D1) .LE. 0.50-10)    GOTO    1 

2  D1=D1*(DA+DFL0AT(K) ) 

4  IF(KOOK.GT.KOUNT)  GOTO  1 
DO  3  L=KOCK.KOUNT 
IF(CDABS(D1) .LE. 0.50-10)  GOTO  1 

3  01  =  D1*(DA  +  DFL0AT(U  ) 
1    Dil=DlH-Dl 

D2=D2-(2.D0/DN) 
DG=OQ*DN*DN 

IFCDQ.    GE. 0.1040)    GOTO    56 

ZZ=(D11    +D1*(DA*DFL0AT(K0UNT) )*D2)*DX/DQ 
KCUNT=KOUNTd 
IFIK0UNT.GT.25)     GOTO    56 
DL=DL+ZZ 

IF(CDABS(ZZ)  .LE. 0.50-8)  GOTO  21 
DN=DN+  l.DO 
GOTO  2  2 

21  RETURN 

56  WRITE(6,53)  11 
GOTO  21 

51  FGRMATdH  ,4 ( 5X , E 15 . 8  )  ) 

53  FGRMATdH  ,35HG(X)  DID  NOT  CONVERGE .  LAST  TERM  WAS,E15. 

END 

SUBRCUTINE  DRDL AP { A, X , RESULT ) 
C0MPLEX*16  A, X, RESULT ,YOA,YDD 
REAL*8  YDB,YDN 

DRDLAP  SOLVES  THE  FIRST  DERIVATIVE  OF  THE  LAGUERRE  FUNCTION 
OF  THE  FIRST  KIND. 

YDA=A+1.DC 
RESULT=A 
YDN=1.D0 
YDD=A 

22  YDB=YDN+1 .CO 
YCD=YDD*YDA*X/  (YCN*YDB) 
RESULT=RESULT+YDD 

IF(CDABS(YDD).LT.0.5D-8)  GOTO  21 
YDA=YDA+1.D0 

YCN=YDN+1 .DO 
GLTO  22 
21  RETURN 

23  FORMATdH  ,  ?.  (  5X  ,  E15  .8  )  ) 

54  FGRMAT  ( 1H  ,6020.7) 

END 

SUBROUTINE  OGDLAP ( A, X , RESULT, PLG, RESUL ) 

COMPLEX* 16  A,X,RESULT,PLG,RESUL,XA,X2,XX,XZ,X4,X3 

REAL*8  XN 

DGDLAP  SCLVES  THE  FIRST  DERIVATIVE  OF  THE  SECOND  KIND  OF 
LAGUERRE  FUNCT ION. WHEN  X=0,  THIS  FUNCTION  DOES  NOT  EXIST. 

IT (CDABS(X).EQ.O.CO)     GOTO    10 

XN=2.D0 

XA=A 

RESUL=RESULT*CDLOG(X)+PLG/Xd .DO- < 2 .D0*X A ) 

X2=1.D0/XA-2.D0 
XX  =  XA 
21    XA=XA+1.D0 

IF(CDABS(XA) .LE.0.5D-4I    GOTO    22 
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XX=XX*X/(XN*XN~XN)*XA 

X2=X2  -2.DC/XN  +1.D0/XA 

XZ=XX*X2 

RESUl=RCSUL  +  XZ 

IF(CDABS(XZ) .LE.0.5D-8)  GOTO  10 

XN=XN  +  1.D0 

IF(CDABS(XZ) .GF.0.1D60)  GOTO  23 

GO  TO  21 

22  KISS=1 

XN=2.D0 
XA=A 

RESUL=RFSULT*CDLOG(X)+PLG/X+l.DO-(2.D0*XA) 
XQ=1.D0 
XX=1.D0 
X2=-2.D0 

5  XX=XX*X 
X3=0.D0 

DO  6  II=OfKISS 

X4=1.D0 

KIN=KISS+1 

DC  7  JJ=0,  II 

IF( JJ.EQ.II)  GOTO  9 

IF(CDABS(X4) .LE.0.5D-10)  GOTO  6 

7  X4=X4*(XA+DFL0AT( JJ) > 

9  IF(KIN.GT.KISS)  GOTO  6 
DC  8  LL=KIN,KISS 
IFtCDABS(XA) .LE. 0.50-101  GOTO  6 

8  X4=X4*(XA+DFL0AT( LL) ) 

6  X3=X3+X4 
X2=X2-(2.DC/XN) 
XC=XQ*XN*DFLOAT(KISS) 

XZ=(X3+X4*(XA+DFL0AT(KISS) )*X2)*XX/XQ 
KISS=KISS+1 

RESUL=RESUL+XZ 

IF(CDABS(XZ) .LE. 0.50-8)  GOTO  10 

XN=XN+1.D0 

IF(KISS.LT.25)  GOTO  5 

WRITE<6,54)  XZ 

23  N=XN 
WPITE(6,24)  N 

10  RETURN 

24  F0RMAT(40H  OVERFLOW  ABOUT  TO  OCCUR  IN  DGDLAP  AFTER,!  3 
+  7H  TERMS.) 

51  FORMATdH  ,4(  5X  ,  E  15.  8  )  )  i 

54  FURMATUH  ,36HG'lX)  DID  NOT  CONVERGE .  LAST  TERM  WAS, 
,  E15.8) 

END 
//GC.FT06F001  DC  DC E= ( R ECFM  =  FA , BLKS I ZE=  1  33) , SP ACE= ( CYL , ( 15, 1 
//GO.SYSUDUMP  DD  SYSOUT=A 
//GO.SYSIN  DD  * 
0.100D  08  0.520D  07  0.729D-04  0. 2000-02  0.500D-01 
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